## Figure 2: RDD analysis

library(dplyr)
library(ggplot2)
library(reshape2)
library(grid)
library(rdlocrand)
library(lubridate)

## Clear workspace

rm(list=ls())

##Opening the voting data:

load("ReplicationData.RData")


################
#### 2008 ######
################

## Turnout variable
data2008$turnout <- as.numeric(data2008$VOTOU.1º.TURNO) - 1


#############################
## Discontinuity analysis ###
#############################

## 18-year olds

## Election Day: subsetting the 14-day window

data.18.rddED <- dplyr::filter(data2008, (dob >= "1990-09-29" & dob <= "1990-10-12") & turnout == 1)

DataRDDED18 <- data.18.rddED %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDED18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDED18$WeekendHolidays <- as.numeric(weekdays(DataRDDED18$dob) %in% c("Saturday", "Sunday") |
                                          DataRDDED18$dob=="1990-10-12")

##Excluding weekend days:

DataRDDED18nowe <- DataRDDED18[ ! DataRDDED18$WeekendHolidays==1,]

##Using the rdlocrand package

rddED18nowe <- rdrandinf(DataRDDED18nowe$voters, DataRDDED18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateED2008 <- rddED18nowe$obs.stat/rddED18nowe$sumstats[3,1]*100
estimateED2008

## % increase: lower bound
estimateED2008lw <- rddED18nowe$ci[1]/rddED18nowe$sumstats[3,1]*100
estimateED2008lw

## % increase: upper bound
estimateED2008up <- rddED18nowe$ci[2]/rddED18nowe$sumstats[3,1]*100
estimateED2008up


## End-of-Year: subsetting the 14-day window

data.18.rddEoY <- dplyr::filter(data2008, (dob >= "1990-12-25" & dob <= "1991-01-07") & turnout == 1)

DataRDDEoY18 <- data.18.rddEoY %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDEoY18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDEoY18$WeekendHolidays <- as.numeric(weekdays(DataRDDEoY18$dob) %in% c("Saturday", "Sunday") | DataRDDEoY18$dob=="1990-12-24" | DataRDDEoY18$dob=="1990-12-25" | DataRDDEoY18$dob=="1990-12-31" | DataRDDEoY18$dob=="1991-01-01")

##Excluding weekend days:

DataRDDEoY18nowe <- DataRDDEoY18[ ! DataRDDEoY18$WeekendHolidays==1,]

##Using the rdlocrand package

rddEoY18nowe <- rdrandinf(DataRDDEoY18nowe$voters, DataRDDEoY18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateEoY2008 <- rddEoY18nowe$obs.stat/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2008

## % increase: lower bound
estimateEoY2008lw <- rddEoY18nowe$ci[1]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2008lw

## % increase: upper bound
estimateEoY2008up <- rddEoY18nowe$ci[2]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2008up


################
#### 2010 ######
################

## Turnout variable
data2010$turnout <- as.numeric(data2010$VOTOU.1º.TURNO) - 1


#############################
## Discontinuity analysis ###
#############################

## 18-year olds

## Election Day: subsetting the 14-day window

data.18.rddED <- dplyr::filter(data2010, (dob >= "1992-09-27" & dob <= "1992-10-10") & turnout == 1)

DataRDDED18 <- data.18.rddED %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDED18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDED18$WeekendHolidays <- as.numeric(weekdays(DataRDDED18$dob) %in% c("Saturday", "Sunday"))

##Excluding weekend days:

DataRDDED18nowe <- DataRDDED18[ ! DataRDDED18$WeekendHolidays==1,]

##Using the rdlocrand package

rddED18nowe <- rdrandinf(DataRDDED18nowe$voters, DataRDDED18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateED2010 <- rddED18nowe$obs.stat/rddED18nowe$sumstats[3,1]*100
estimateED2010

## % increase: lower bound
estimateED2010lw <- rddED18nowe$ci[1]/rddED18nowe$sumstats[3,1]*100
estimateED2010lw

## % increase: upper bound
estimateED2010up <- rddED18nowe$ci[2]/rddED18nowe$sumstats[3,1]*100
estimateED2010up


## End-of-Year: subsetting the 14-day window

data.18.rddEoY <- dplyr::filter(data2010, (dob >= "1992-12-25" & dob <= "1993-01-07") & turnout == 1)

DataRDDEoY18 <- data.18.rddEoY %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDEoY18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDEoY18$WeekendHolidays <- as.numeric(weekdays(DataRDDEoY18$dob) %in% c("Saturday", "Sunday") | DataRDDEoY18$dob=="1992-12-24" | DataRDDEoY18$dob=="1992-12-25" | DataRDDEoY18$dob=="1992-12-31" | DataRDDEoY18$dob=="1993-01-01")

##Excluding weekend days:

DataRDDEoY18nowe <- DataRDDEoY18[ ! DataRDDEoY18$WeekendHolidays==1,]

##Using the rdlocrand package

rddEoY18nowe <- rdrandinf(DataRDDEoY18nowe$voters, DataRDDEoY18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateEoY2010 <- rddEoY18nowe$obs.stat/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2010

## % increase: lower bound
estimateEoY2010lw <- rddEoY18nowe$ci[1]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2010lw

## % increase: upper bound
estimateEoY2010up <- rddEoY18nowe$ci[2]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2010up


################
#### 2012 ######
################

## Turnout variable
data2012$turnout <- as.numeric(data2012$VOTOU.1º.TURNO) - 1


## 18-year olds

## Election Day: subsetting the 14-day window

data.18.rddED <- dplyr::filter(data2012, (dob >= "1994-10-01" & dob <= "1994-10-14") & turnout == 1)

DataRDDED18 <- data.18.rddED %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDED18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDED18$WeekendHolidays <- as.numeric(weekdays(DataRDDED18$dob) %in% c("Saturday", "Sunday") | DataRDDED18$dob=="1994-10-12")

##Excluding weekend days:

DataRDDED18nowe <- DataRDDED18[ ! DataRDDED18$WeekendHolidays==1,]

##Using the rdlocrand package

rddED18nowe <- rdrandinf(DataRDDED18nowe$voters, DataRDDED18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateED2012 <- rddED18nowe$obs.stat/rddED18nowe$sumstats[3,1]*100
estimateED2012

## % increase: lower bound
estimateED2012lw <- rddED18nowe$ci[1]/rddED18nowe$sumstats[3,1]*100
estimateED2012lw

## % increase: upper bound
estimateED2012up <- rddED18nowe$ci[2]/rddED18nowe$sumstats[3,1]*100
estimateED2012up


## End-of-Year: subsetting the 14-day window

data.18.rddEoY <- dplyr::filter(data2012, (dob >= "1994-12-25" & dob <= "1995-01-07") & turnout == 1)

DataRDDEoY18 <- data.18.rddEoY %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDEoY18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDEoY18$WeekendHolidays <- as.numeric(weekdays(DataRDDEoY18$dob) %in% c("Saturday", "Sunday") | DataRDDEoY18$dob=="1994-12-24" | DataRDDEoY18$dob=="1994-12-25" | DataRDDEoY18$dob=="1994-12-31" | DataRDDEoY18$dob=="1995-01-01")

##Excluding weekend days:

DataRDDEoY18nowe <- DataRDDEoY18[ ! DataRDDEoY18$WeekendHolidays==1,]

##Using the rdlocrand package

rddEoY18nowe <- rdrandinf(DataRDDEoY18nowe$voters, DataRDDEoY18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateEoY2012 <- rddEoY18nowe$obs.stat/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2012

## % increase: lower bound
estimateEoY2012lw <- rddEoY18nowe$ci[1]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2012lw

## % increase: upper bound
estimateEoY2012up <- rddEoY18nowe$ci[2]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2012up


################
#### 2014 ######
################

## Turnout variable
data2014$turnout <- as.numeric(data2014$VOTOU.1º.TURNO) - 1


#############################
## Discontinuity analysis ###
#############################

## 18-year olds

## Election Day: subsetting the 14-day window

data.18.rddED <- dplyr::filter(data2014, (dob >= "1996-09-29" & dob <= "1996-10-12") & turnout == 1)

DataRDDED18 <- data.18.rddED %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDED18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDED18$WeekendHolidays <- as.numeric(weekdays(DataRDDED18$dob) %in% c("Saturday", "Sunday") |
                                          DataRDDED18$dob=="1996-10-12")

##Excluding weekend days:

DataRDDED18nowe <- DataRDDED18[ ! DataRDDED18$WeekendHolidays==1,]

##Using the rdlocrand package

rddED18nowe <- rdrandinf(DataRDDED18nowe$voters, DataRDDED18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateED2014 <- rddED18nowe$obs.stat/rddED18nowe$sumstats[3,1]*100
estimateED2014

## % increase: lower bound
estimateED2014lw <- rddED18nowe$ci[1]/rddED18nowe$sumstats[3,1]*100
estimateED2014lw

## % increase: upper bound
estimateED2014up <- rddED18nowe$ci[2]/rddED18nowe$sumstats[3,1]*100
estimateED2014up


## End-of-Year: subsetting the 14-day window

data.18.rddEoY <- dplyr::filter(data2014, (dob >= "1996-12-25" & dob <= "1997-01-07") & turnout == 1)

DataRDDEoY18 <- data.18.rddEoY %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDEoY18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDEoY18$WeekendHolidays <- as.numeric(weekdays(DataRDDEoY18$dob) %in% c("Saturday", "Sunday") | DataRDDEoY18$dob=="1996-12-24" | DataRDDEoY18$dob=="1996-12-25" | DataRDDEoY18$dob=="1996-12-31" | DataRDDEoY18$dob=="1997-01-01")

##Excluding weekend days:

DataRDDEoY18nowe <- DataRDDEoY18[ ! DataRDDEoY18$WeekendHolidays==1,]

##Using the rdlocrand package

rddEoY18nowe <- rdrandinf(DataRDDEoY18nowe$voters, DataRDDEoY18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateEoY2014 <- rddEoY18nowe$obs.stat/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2014

## % increase: lower bound
estimateEoY2014lw <- rddEoY18nowe$ci[1]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2014lw

## % increase: upper bound
estimateEoY2014up <- rddEoY18nowe$ci[2]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2014up


################
#### 2016 ######
################

## Turnout variable
data2016$turnout <- as.numeric(data2016$VOTOU.1º.TURNO) - 1


#############################
## Discontinuity analysis ###
#############################

## 18-year olds

## Election Day: subsetting the 14-day window

data.18.rddED <- dplyr::filter(data2016, (dob >= "1998-09-26" & dob <= "1998-10-09") & turnout == 1)

DataRDDED18 <- data.18.rddED %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDED18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDED18$WeekendHolidays <- as.numeric(weekdays(DataRDDED18$dob) %in% c("Saturday", "Sunday"))

##Excluding weekend days:

DataRDDED18nowe <- DataRDDED18[ ! DataRDDED18$WeekendHolidays==1,]

##Using the rdlocrand package

rddED18nowe <- rdrandinf(DataRDDED18nowe$voters, DataRDDED18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateED2016 <- rddED18nowe$obs.stat/rddED18nowe$sumstats[3,1]*100
estimateED2016

## % increase: lower bound
estimateED2016lw <- rddED18nowe$ci[1]/rddED18nowe$sumstats[3,1]*100
estimateED2016lw

## % increase: upper bound
estimateED2016up <- rddED18nowe$ci[2]/rddED18nowe$sumstats[3,1]*100
estimateED2016up


## End-of-Year: subsetting the 14-day window

data.18.rddEoY <- dplyr::filter(data2016, (dob >= "1998-12-25" & dob <= "1999-01-07") & turnout == 1)

DataRDDEoY18 <- data.18.rddEoY %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDEoY18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDEoY18$WeekendHolidays <- as.numeric(weekdays(DataRDDEoY18$dob) %in% c("Saturday", "Sunday") | DataRDDEoY18$dob=="1998-12-24" | DataRDDEoY18$dob=="1998-12-25" | DataRDDEoY18$dob=="1998-12-31" | DataRDDEoY18$dob=="1999-01-01")

##Excluding weekend days:

DataRDDEoY18nowe <- DataRDDEoY18[ ! DataRDDEoY18$WeekendHolidays==1,]

##Using the rdlocrand package

rddEoY18nowe <- rdrandinf(DataRDDEoY18nowe$voters, DataRDDEoY18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateEoY2016 <- rddEoY18nowe$obs.stat/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2016

## % increase: lower bound
estimateEoY2016lw <- rddEoY18nowe$ci[1]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2016lw

## % increase: upper bound
estimateEoY2016up <- rddEoY18nowe$ci[2]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2016up


################
#### 2018 ######
################

## Turnout variable
data2018$turnout <- as.numeric(data2018$VOTOU.1º.TURNO) - 1


#############################
## Discontinuity analysis ###
#############################

## 18-year olds

## Election Day: subsetting the 14-day window

data.18.rddED <- dplyr::filter(data2018, (dob >= "2000-10-01" & dob <= "2000-10-14") & turnout == 1)

DataRDDED18 <- data.18.rddED %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDED18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDED18$WeekendHolidays <- as.numeric(weekdays(DataRDDED18$dob) %in% c("Saturday", "Sunday") |
                                          DataRDDED18$dob=="2000-10-12")

##Excluding weekend days:

DataRDDED18nowe <- DataRDDED18[ ! DataRDDED18$WeekendHolidays==1,]

##Using the rdlocrand package

rddED18nowe <- rdrandinf(DataRDDED18nowe$voters, DataRDDED18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateED2018 <- rddED18nowe$obs.stat/rddED18nowe$sumstats[3,1]*100
estimateED2018

## % increase: lower bound
estimateED2018lw <- rddED18nowe$ci[1]/rddED18nowe$sumstats[3,1]*100
estimateED2018lw

## % increase: upper bound
estimateED2018up <- rddED18nowe$ci[2]/rddED18nowe$sumstats[3,1]*100
estimateED2018up


## End-of-Year: subsetting the 14-day window

data.18.rddEoY <- dplyr::filter(data2018, (dob >= "2000-12-25" & dob <= "2001-01-07") & turnout == 1)

DataRDDEoY18 <- data.18.rddEoY %>%
  filter(GRAU.INSTRUÇÃO != "Analfabeto") %>%
  group_by(dob) %>%
  summarise(voters = n(),
            turnout = unique(turnout))

DataRDDEoY18$daysToFrom <- seq(6, -7, by = -1)

##Identifying weekend days and holidays:

DataRDDEoY18$WeekendHolidays <- as.numeric(weekdays(DataRDDEoY18$dob) %in% c("Saturday", "Sunday") | DataRDDEoY18$dob=="2000-12-24" | DataRDDEoY18$dob=="2000-12-25" | DataRDDEoY18$dob=="2000-12-31" | DataRDDEoY18$dob=="2001-01-01")

##Excluding weekend days:

DataRDDEoY18nowe <- DataRDDEoY18[ ! DataRDDEoY18$WeekendHolidays==1,]

##Using the rdlocrand package

rddEoY18nowe <- rdrandinf(DataRDDEoY18nowe$voters, DataRDDEoY18nowe$daysToFrom, wl = -7, wr = 6, seed = 50, ci = .05)

## % increase
estimateEoY2018 <- rddEoY18nowe$obs.stat/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2018

## % increase: lower bound
estimateEoY2018lw <- rddEoY18nowe$ci[1]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2018lw

## % increase: upper bound
estimateEoY2018up <- rddEoY18nowe$ci[2]/rddEoY18nowe$sumstats[3,1]*100
estimateEoY2018up



## Figure 2: plotting the effects with CI

setEPS()
postscript("Figure_2_Effects_of_CV.eps",  width = 7.5, height = 4.5)

par(mfrow = c(1, 2), mar = c(1, 4, 1, 1), las = 1)

plot(c(0.1, 0.2, 0.3), c(estimateED2010, estimateED2014, estimateED2018), pch = c(16), xlim = c(0, 0.8), ylim = c(-4.2, 48), 
     xlab = "", ylab = "Estimate increase in %", main = "", axes = FALSE)

points(c(0.5, 0.6, 0.7), c(estimateEoY2010, estimateEoY2014, estimateEoY2018), pch = c(17))

axis(side = 2, at = seq(from = -2, to = 42, by = 4), cex.axis=0.79, las=1)

text(.1, -4.2, "2010", cex = .6, font = 2)
text(.2, -4.2, "2014", cex = .6, font = 2)
text(.3, -4.2, "2018", cex = .6, font = 2)
text(.5, -4.2, "2010", cex = .6, font = 2)
text(.6, -4.2, "2014", cex = .6, font = 2)
text(.7, -4.2, "2018", cex = .6, font = 2)

lines(c(0.1, 0.1), c(estimateED2010lw, estimateED2010up)) 
lines(c(0.2, 0.2), c(estimateED2014lw, estimateED2014up)) 
lines(c(0.3, 0.3), c(estimateED2018lw, estimateED2018up)) 

lines(c(0.5, 0.5), c(estimateEoY2010lw, estimateEoY2010up)) 
lines(c(0.6, 0.6), c(estimateEoY2014lw, estimateEoY2014up)) 
lines(c(0.7, 0.7), c(estimateEoY2018lw, estimateEoY2018up)) 

abline(h = 0, lty = "dashed")

text(.2, 44, "Election Day", cex = 1, font = 1)
text(.6, 44, "End-of-year", cex = 1, font = 1)

text(.4, 47, "Presidential elections", cex = 1.2, font = 2)


plot(c(0.1, 0.2, 0.3), c(estimateED2008, estimateED2012, estimateED2016), pch = c(16), xlim = c(0, 0.8), ylim = c(-4.2, 48), 
     xlab = "", ylab = "Estimate increase in %", main = "", axes = FALSE)

points(c(0.5, 0.6, 0.7), c(estimateEoY2008, estimateEoY2012, estimateEoY2016), pch = c(17))

axis(side = 2, at = seq(from = -2, to = 42, by = 4), cex.axis=0.79, las=1)

text(.1, -4.2, "2008", cex = .6, font = 2)
text(.2, -4.2, "2012", cex = .6, font = 2)
text(.3, -4.2, "2016", cex = .6, font = 2)
text(.5, -4.2, "2008", cex = .6, font = 2)
text(.6, -4.2, "2012", cex = .6, font = 2)
text(.7, -4.2, "2016", cex = .6, font = 2)

lines(c(0.1, 0.1), c(estimateED2008lw, estimateED2008up)) 
lines(c(0.2, 0.2), c(estimateED2012lw, estimateED2012up)) 
lines(c(0.3, 0.3), c(estimateED2016lw, estimateED2016up)) 

lines(c(0.5, 0.5), c(estimateEoY2008lw, estimateEoY2008up)) 
lines(c(0.6, 0.6), c(estimateEoY2012lw, estimateEoY2012up)) 
lines(c(0.7, 0.7), c(estimateEoY2016lw, estimateEoY2016up)) 

abline(h = 0, lty = "dashed")

text(.2, 44, "Election Day", cex = 1, font = 1)
text(.6, 44, "End-of-year", cex = 1, font = 1)

text(.4, 47, "Municipal elections", cex = 1.2, font = 2)

dev.off()


